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Overview 

Nanothermodynamics  Applied  to  Thermal  Processes  in  Heterogeneous  Materials 

Small- system  thermodynamics  is  used  as  a  guide  to  investigate  some  fundamental 
issues  in  the  behavior  of  complex  systems.  This  “nanothermodynamics”  provides  a  novel 
paradigm  for  interpreting  the  nonlinear  and  inhomogeneous  response  inside  bulk 
materials,  leading  to  improved  modeling  of  the  response.  Results  from  computer 
simulations  of  standard  models  agree  favorably  with  a  wide  range  of  measurements  of 
local  thermal  and  dynamic  properties.  Progress  in  understanding  basic  thermodynamic 
principles  will  enhance  the  predictive  power  and  accuracy  of  most  computer  simulations. 

Unlike  the  usual  theoretical  limit  of  linear  and  homogeneous  behavior  for  large 
systems,  the  laws  of  thermodynamics  are  more  restrictive  when  applied  to  small  systems. 
For  example,  local  entropic  forces  often  yield  non-extensive  changes  in  energy  during 
normal  thermal  fluctuations,  which  must  be  added  to  the  interaction  if  total  energy  is  to 
be  conserved.  Similarly,  local  configurational  entropies  often  yield  nonlinear  terms  in  the 
reaction  rate,  which  must  be  included  in  the  dynamics  if  total  entropy  is  to  be  maximized. 
Nanothermodynamics  helps  to  identify  various  nonlinear  corrections  that  occur  in  the 
dynamics  of  small  systems. 

Finite-size  effects  yield  nonlinear  corrections  due  to  changes  in  the  local  entropy 
during  normal  thermal  fluctuations.  Fluctuations  in  the  configurational  entropy  come 
from  changes  in  the  local  alignment  of  dipoles,  or  the  local  density  of  particles,  which 
can  alter  reaction  rates  by  several  orders  of  magnitude  across  length  scales  of  nanometers. 
Using  computer  simulations  of  standard  models,  such  as  the  Ising  ferromagnet,  the 
nonlinear  corrections  are  shown  to  restore  conservation  of  energy  and  maintain  maximum 
entropy.  Furthermore,  the  nonlinear  correction  improves  agreement  between  computer 
simulations  and  the  measured  properties  of  many  materials.  A  general  result  is  that  finite- 
size  effects  cause  simple  models  to  exhibit  complex  thermal  and  dynamic  behavior. 

Useful  implementation  of  our  basic  research  requires  adapting  the  fundamental 
principles  we  investigate  to  help  guide  advanced  applications.  We  now  collaborate  with 
experts  in  the  area  of  energetic  materials  at  the  University  of  Missouri.  One  result  is  an 
analysis  of  finite-size  thermal  effects  in  molecular-dynamics  simulations  of  nitromethane. 
Future  applications  of  nanothermodynamics  will  be  facilitated  by  continuing  to  exchange 
knowledge  with  such  experts.  Our  goal  is  to  optimize  the  accuracy  and  efficiency  of 
large-scale  simulations  of  the  complex  dynamics  in  energetic  materials. 
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I.  Introduction 

Computer  simulations  are  often  the  best  way  to  investigate  microscopic  details  in 
the  complex  behavior  of  energetic  materials. 1  For  example,  classical  molecular  dynamics 
(MD)  simulations  allow  large-scale  and  long-time  studies  of  atomic-level  motion, 
revealing  details  about  the  structural  and  chemical  consequences  of  shock  waves,"  and 
the  dynamics  of  melting/1  We  explore  finite-size  effects  in  thermodynamics.  Specifically, 
we  investigate  nonlinear  corrections  to  Boltzmann  statistics  that  come  from  local  entropic 
forces.  Our  goal  is  to  fully  understand  the  statistical  thermodynamics  of  finite-sized 
systems,  which  will  enhance  the  efficiency  and  accuracy  of  computer  simulations. 

MD  simulations  of  large-scale  models  give  good  agreement  with  several  features 
in  the  response  of  energetic  materials,  including  the  density,  pressure,  velocity,  and 
temperature  of  detonation  waves.4  However,  most  computer  simulations  exhibit  non- 
extensive  energy  fluctuations  in  both  MD5  and  Monte-Carlo6  algorithms.  In  our  research 
we  emphasize  that  these  simulations  violate  the  basic  laws  of  thermodynamics,  and  we 
investigate  ways  to  improve  the  algorithms.  First  I  will  present  evidence  that  a  more- 
sophisticated  model  that  is  commonly  used  to  simulate  energetic  materials  also  exhibits 
some  fundamental  problems  in  its  thermodynamic  behavior,  then  I  will  summarize  our 
efforts  to  use  nanothermodynamics  as  a  guide  to  address  these  issues. 

II.  Fundamental  issues  in  the  thermal  behavior  of  molecular-dynamics  simulations 

We  analyzed  size-dependent  energy  fluctuations  from  molecular-dynamics 
simulations  of  liquid  nitromethane.  The  system  was  simulated  at  a  temperature  of  360  K 
using  the  Sandia  National  Laboratory  program  LAMMPS  (Large-scale  Atomic/Molecular 
Massively  Parallel  Simulator).  The  simulations  were  kindly  provided  by  Dr.  Lan  He  and 
Prof.  Tommy  Sewell  from  the  University  of  Missouri.  We  used  three  files  from  each 
simulation.  The  files  give  the  position,  velocity,  and  force  on  each  atom  at  a  series  of 
1200  time  steps.  The  data  were  converted  into  the  kinetic  energy  (K),  potential  energy 
change  (AU),  and  total  energy  (E=K+AU)  for  each  atom  at  each  time  step.  These 
simulations  used  1896  nitromethane  molecules  (13272  atoms)  in  a  total  volume  of: 
(2*27.9126  A)  .  The  size-dependent  behavior  was  analyzed  by  subdividing  the  total 
volume  into  2  ,  3  ,  4 , 6  ,  8',  or  12  cubic-shaped  regions;  yielding  sides  for  each  region 
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ranging  in  length  from  L=27. 9126/6  A  to  27.9126  A.  With  a  total  number  of  N(L)  atoms 
in  each  region,  the  values  obtained  were  the  kinetic  energy  K(L),  potential  energy  change 
AU(L),  and  total  energy  E(L).  Specifically,  the  kinetic  energy  was  found  using  K(L)=  Vi 
£j  m,  (vix  +  Viy"+ViZ ),  where  m,  is  the  atomic  mass  of  each  atom,  and  the  sum  is  over  all 
atoms  in  the  region.  A  factor  of  2390.06  (kcal/g)/(A/fs)  was  used  to  convert  the 
simulation  units  (g/mole)*(A/fs)  into  kcal/mole.  The  potential  energy  change  was 
calculated  by  integrating  the  power  over  each  time  step,  U(L)=/(F*v)dt,  using  a  simple 
trapezoid  integration.  A  factor  of  75  fs/(time-step)  was  used  to  convert  the  simulation 
units  [kcal/(mole-A)]*(A/fs)  into  kcal/mole.  The  energies  were  averaged  over  the  time 
steps,  and  divided  by  <N(L)>  to  obtain  the  per  mole  values:  <k(L)>,  <u(L)>,  and  <z{L)>. 

In  a  similar  manner  the  average  fluctuations  per  particle  were  also  obtained:  <(Ak(L))“>, 

2  2 

<(Au(L))->,  and  <(Ae(L))“>.  Figures  1  and  2  show  some  results. 

Energy  fluctuations  yield  fundamental  thermodynamic  quantities  of  heat  capacity 

and  entropy.  Specifically,  cy  =  <((AE/T)2>/(knN)  gives  the  specific  heat  in  units  of 

Boltzmann’s  constant  kB,  while  A sl  =  \(cvlT)dT  gives  the  change  in  entropy  per  particle 

between  two  temperatures.  The  equipartition  theorem  states  that  all  classical  particles 

have  a  heat  capacity  of  Vika  per  degree  of  freedom;  but  the  fluctuations  shown  in  Fig.  2 

decrease  with  decreasing  system  size.  Furthermore,  small-system  thermodynamics 

(nano thermodynamics)  requires  that  the  entropy  is  additive  and  extensive,  unlike  the 

nonextensive  entropy  deduced  from  Fig.  2.  Thus,  classical  MD  simulations  of  energetic 

materials  may  exhibit  unphysical  behavior  on  small  length  scales. 

Figure  1  Fluctuations  in  kinetic 
energy  per  mole  as  a  function  of  the 
side-length  (in  Angstroms)  of  cube¬ 
shaped  regions.  The  results  are 
from  MD  simulations  of  liquid 
nitromethane  at  a  temperature  of 
360K.  Symbols  and  error  bars  are 
from  the  average  and  standard 
deviation  over  the  interval  from 
600-1200  time  steps.  The  relatively 
large  error  bars  prohibit  a  clear 
choice  for  the  functional  form  of 
the  size  dependence. 

2  4  6  8  10  12 

L  (A) 


_Q) 

O 


03  250 

o 

A 


200 

<1 

V 

i  Rn 


n  1 i  1 i  1 i  1 i  1 r 

i  MD  simulation  Nitromethane  360  K 


.  *  t 


I 


_l _ i _ I _ i _ I _ i _ I _ i _ I _ i _ L_ 


2 


Figure  2  Fluctuations  in  potential 
energy  per  mole  as  a  function  of 
the  length  (in  Angstroms)  of  cube¬ 
shaped  internal  regions.  Symbols 
and  error  bars  are  from  the 
average  and  standard  deviation 
over  the  interval  from  1000-1200 
time  steps.  The  solid  curve  is  a  fit 
to  the  data  showing  that  the 
fluctuations  decrease  inversely 
proportional  to  the  length  of  the 
region. 

fluctuations  in  MD  simulations  of 

simple  fluids  using  the  Lennard-Jones  (L-J)  model,  as  shown  in  Fig.  3.  Specifically,  for 

fluctuations  in  the  potential  energy  at  low  temperatures  where  the  system  is  liquid- like 

(red  circles),  the  resulting  specific  heat  approaches  the  theoretically  expected  value  of 

cv/kB~> 3/2  for  large  regions,  but  cy/kB  decreases  inversely  proportional  to  the  region 

length  (red-dashed  curve).  Fluctuations  in  kinetic  energy  (black  squares)  also  yield  cv/kB 

3/2  for  large  regions,  but  this  cylkB  decreases  inversely  proportional  to  volume  for 

small  regions  (black-solid  curve).  Thus  the  L-J  model  provides  a  simplified  system  to 

study  basic  issues  that  occur  in  most  MD  simulations.  The  goal  of  our  current  research  is 

to  understand  why  these  MD  simulations  show  fundamentally  non-physical  behavior,  and 

how  can  they  be  corrected  to  exhibit  more  realistic  behavior,  consistent  with  the  laws  of 

Figure  3  Specific  heat  as  a  function  of  length  (in 
units  of  the  L-J  lattice)  of  cube-shaped  regions 
from  MD  simulations  of  the  L-J  fluid.  Colored 
symbols  are  from  fluctuations  in  potential 
energy  per  particle  at  five  reduced  temperatures 
(given  in  legend).  Black  squares  are  from 
fluctuations  in  kinetic  energy,  after  averaging 
over  all  five  temperatures  (error  bars  are 
comparable  to  symbol  size).  The  lines  give  a 
best  fit  to  the  data  using  an  inverse  volume 
dependence  for  the  kinetic  energy  (solid  black), 
and  an  inverse  length  dependence  for  the 
potential  energy  (dashed  colored).  Note  that  both 
the  solid-black  and  dashed-red  lines  yield  Cy/kB 
—3/2  for  large  regions  (values  are  given  in  the 
legend),  as  expected  by  the  equipartition 
theorem.  However,  the  simulations  all  show  a 
decrease  in  Cy/kB  for  small  regions,  similar  to 
that  shown  in  Fig.  2  for  nitromethane.  This  size- 
dependent  cv/kB  violates  the  equipartition 
theorem,  and  the  thermodynamic  requirement 
that  entropy  is  extensive. 
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III.  Fundamental  issues  in  the  thermal  behavior  of  Monte- Carlo  simulations 


A  major  achievement  during  the  just-completed  period  of  ARO  support  was  to 
characterize  and  correct  violations  of  fundamental  thermodynamic  behavior  in  Monte- 
Carlo  (MC)  simulations  of  the  Ising  model.7  The  solid  black  lines  in  Fig.  4  show  results 
using  the  uncorrected  (Metropolis)  algorithm  on  the  standard  Ising  model.  The  energy  per 
particle  in  Fig.  4  a)  is  independent  of  system  size,  whereas  the  energy  fluctuations  per 
particle  in  Fig.  4  b)  decrease  with  decreasing  region  size,  similar  to  that  shown  in  Figs.  2 
and  3  for  MD  simulations.  Again,  because  the  energy  fluctuations  yield  the  heat  capacity 
and  change  in  entropy  between  temperatures,  the  behavior  in  Fig.  4  b)  violates  the 
thermodynamic  principle  that  entropy  should  be  extensive. 

The  Metropolis  algorithm  that  is  used  for  most  MC  simulations  is  based  on  the 
Boltzmann  factor  e~AUIkBT  .  Although  this  is  the  most  widely  used  formula  in  all  of 
statistical  mechanics,  it  is  based  on  some  restrictive  assumptions:  the  system  must  couple 

o 

instantly  to  a  large  heat  bath,  and  heat  must  be  the  only  source  of  entropy.  In  the  past  20 
years,  several  experimental  techniques  have  shown  that  the  primary  response  of  most 
materials  comes  from  an  ensemble  of  effectively  independent  nanometer- sized  regions. 
Thus,  local  thermal  properties  play  an  important  role  in  the  primary  response.  Indeed,  the 
Boltzmann  factor  is  strictly  valid  only  for  independent  particles  that  couple  individually 


to  an  infinite  heat  bath,  which  is  never  true  for  classical  systems  of  interacting  particles. 
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Figure  4  (a)  Average  interaction  energy  per 
particle  and  (b)  specific  heat  as  a  function  of 
region  size  at  two  temperatures,  from  ref.  7. 
Simulations  are  performed  on  an  Ising  system  of 
N= 27,000  particles  with  periodic  boundary 
conditions.  The  system  is  subdivided  into  140 
regions,  ranging  in  size  from  n= 6  to  1080 
particles.  The  solid  (black)  lines  come  from  the 
standard  Metropolis  algorithm  with  no 
correction,  i.e.  g= 0  where  g  is  the  nonlinear 
constraint  parameter.  Note  that  these  uncorrected 
simulations  show  a  reduction  in  energy 
fluctuations  for  small  regions,  similar  to  the 
behavior  shown  in  Figs.  2  and  3  for  molecular 
dynamics  simulations.  Symbols  come  from 
simulations  with  nonlinear  correction  values 
given  in  the  legend.  The  dashed  lines  show  that 
for  g=  1 ,  in  (a)  the  energy  per  particle  is  reduced 
inversely  proportional  to  region  size  <Uln>~\ln, 
and  in  (b)  the  specific  heat  is  independent  of 
region  size.  Error  bars  are  shown  for  the  g=l 
data. 
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IV.  Beyond  the  Boltzmann  factor 

We  improve  the  thermal  and  dynamic  properties  of  the  Ising  model  by  modifying 
the  Boltzmann  factor,  essentially  providing  a  nonlinear  correction  to  temperature.  The 
correction  comes  from  the  entropy  of  small  regions  in  the  sample,  which  have  finite-size 
effects  that  are  not  in  standard  Boltzmann  statistics.  Indeed,  the  Boltzmann  factor  gives 
the  thermal  activation  rate  only  if  the  local  entropy  is  at  its  peak  value;  but  local  regions 
often  fluctuate  into  a  low-entropy  state  with  dynamics  that  is  significantly  slower. 
Furthermore,  we  have  found  that  these  nonlinear  corrections  are  necessary  for 
conservation  of  energy  in  small  systems,  and  for  the  uniform  specific  heat  that  resolves 

Q 

Gibbs’  paradox  for  interacting  particles.  Moreover,  we  have  shown  that  the  correction 
factor  greatly  improves  the  agreement  between  simulations  of  the  Ising  model  and  the 
measured  response  of  ferromagnetic  materials  and  critical  fluids,9  as  shown  in  Fig.  5.  The 
corrections  are  also  necessary  for  detailed  balance  between  a  superposition  of  equivalent 
states,  allowing  computer  simulations  of  classical  models  to  exhibit  quantum- like 
statistics.  In  fact  the  nonlinear  correction  yields  an  entropic  force,  similar  to  the  entropic 
force  that  favors  low-energy  states  at  low  temperature.  Thus,  we  use  this  nonlinear 
correction  to  improve  standard  models  of  simplified  systems,  with  the  goal  of  transferring 
our  knowledge  to  help  facilitate  more  accurate  and  efficient  simulations  of  sophisticated 
models  of  energetic  materials. 


Figure  5  Effective  scaling  exponent  yeff=  -dlog(j) 
/  dlog(r)  as  a  function  of  reduced  temperature  r  = 

( T-Tc)/Tc ,  where  /  is  the  susceptibility  and  Tc  is 
the  critical  temperature,  from  Ref.  9.  Symbols  are 
from  measured  susceptibilities  of  two  magnetic 
materials  (Gd-«  from  and  Ni-»),  and  one  critical 
fluid  (CO2- A).  Lines  are  from  simulations  of  the 
Ising  model  using:  the  standard  Boltzmann  factor 
g= 0  (red,  dashed);  and  fluctuation  constraint  with 
constraint  parameter  g=0. 5  (green,  dash-dot),  g=  I 
(black,  solid),  and  g= 2  (blue,  dash-dot-dot).  Note 
the  failure  of  the  standard  Ising  model  (g= 0)  to 
match  measured  susceptibilities. 
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V.  Applications  of  nanothermodynamics 

We  study  fundamental  issues  in  the  statistical  thermodynamics  of  small  systems. 
Our  approach  utilizes  nonlinear  corrections  to  classical  Boltzmann  statistics  that  restore 
extensive  entropy,  conservation  of  energy,  and  detailed  balance  between  a  superposition 
of  equivalent,  quantum- like  states,  as  described  in  the  previous  sections.  Because  these 
corrections  yield  a  type  of  entropic  force,  similar  to  the  force  that  causes  heat  flow  from  a 
thermal  gradient,  they  can  be  thought  of  as  a  nonlinear  correction  to  temperature. 

Temperature  (7)  is  the  familiar  variable  that  governs  thermal  fluctuations  and 
reaction  rates  in  materials,  but  this  T  comes  from  just  the  linear  term  in  the  change  of 
entropy  with  changing  energy;  nonlinear  terms  may  contribute  significantly  to  the  net 
behavior.  The  influence  of  temperature  is  usually  found  from  the  Boltzmann  factor 
P  e~AulkBT  ;  which  gives  the  relative  probability  (e.g.  reaction  rate)  for  a  system  to 
change  its  energy  by  A U.  This  Boltzmann  factor  is  the  most  widely  used  formula  in 
statistical  thermodynamics,  providing  the  basis  for  Maxwell- Boltzmann,  Bose-Einstein, 
and  Fermi-Dirac  statistics.  However,  the  Boltzmann  factor  is  a  linear  approximation, 
valid  only  for  large  systems,  in  thermal  equilibrium  at  low  temperatures,  with  negligible 
fluctuations.  Thus,  for  many  systems  at  normal  temperatures  it  is  necessary  to  use  the 
more  fundamental  expression  P  «=  eASIkB  for  the  probability  of  a  local  process,  where  AS 
is  the  offset  in  entropy  from  thermal  equilibrium  due  to  a  local  fluctuation. 

Nonlinear  corrections  to  temperature  are  necessary  to  simulate  several  features  in 
the  thermal  and  dynamic  properties  of  complex  materials.  As  was  shown  in  Fig.  5,  the 
critical  behavior  of  several  samples  is  improved  using  the  nonlinear  correction;  Fig.  6 
shows  another  example.  The  solid  (black)  line  in  the  main  part  of  Fig.  6  comes  from  the 
excess  specific  heat  of  FaMnCF  measured  near  the  Jahn-Teller  structural  distortion  at 
Tc= 735  K.Error!  Bookmark  not  defined'  The  dashed  (red)  line  comes  from  simulations  of  the 
antiferromagnetic  (AF)  Ising  model  using  standard  Boltzmann  statistics  (Metropolis 
algorithm).  The  solid  squares  come  from  the  same  model,  but  with  a  nonlinear  correction 
from  the  local  entropy  of  3x3x3  regions.  The  inset  of  Fig.  6  shows  the  potential  energy 
per  particle  of  this  AF  Ising  model  as  a  function  of  temperature.  Again  the  dashed  (red) 
line  comes  from  standard  simulations  using  the  Boltzmann  factor  alone.  The  symbols 
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come  from  the  same  model,  but  with  nonlinear  corrections  from  cubic  regions  of  various 


sizes.  Large  regions 


Figure  6.  Main  plot  shows  specific  heat  as  a  function  of 
normalized  temperature.  Solid  line  is  from 
measurements  of  LaMnCL  Dashed  line  is  from 
simulations  of  the  AF  fsing  model  using  the  standard 
Boltzmann  factor.  Solid  squares  with  error  bars  are  from 
our  simulations  of  the  same  model,  but  include  a 
nonlinear  correction  to  temperature  from  the  local 
entropy  of  3x3x3  regions,  with  no  other  adjustable 
parameters.  Inset  shows  the  average  energy  per  particle 
as  a  function  of  temperature.  Again  the  dashed  line  is 
from  simulations  using  the  standard  Boltzmann  factor. 
The  symbols  include  a  nonlinear  correction  factor  from 
the  local  entropy  of  regions  of  size  16x16x16-  , 

12xl2xl2-4,  8x8x8- Y,  6x6x6- A,  4x4x4-*,  and  3x3x3- 
■.  Because  of  negligible  correction  for  large  regions, 
several  symbols  lie  underneath  the  dashed  line. 


(16x16x16  &  12x12x12)  have  negligible  influence  on  the  energy  because  they  are 
unlikely  to  fluctuate  into  a  low-entropy  state.  Small  regions  (3x3x3)  have  large 
fluctuations  that  significantly  lower  the  energy  at  all  temperatures,  while  increasing  the 
transition  temperature  by  nearly  50%.  Also  note  the  sharpening  of  the  transition, 
consistent  with  the  measurements  shown  in  the  main  part  of  Fig.  6.  Direct  evidence  for 
the  existence  of  small  correlated  regions  in  LaMnOs  comes  from  neutron  scattering 
data.10  The  red  circles  in  Fig.  7  show  the  correlation  length  deduced  from  the  pair 
distribution  function  measured  in  LaMnOs.  The  black  squares  come  from  our  simulations 
of  the  AF  Ising  model,  showing  similar  short-range  order  above  Tc.  The  simulations  use 
the  measured  value  of  0.40  nm  for  the  average  lattice  parameter,  and  3x3x3  regions 
consistent  with  the  lowest  energy  shown  in  the  inset  of  Fig.  6,  with  no  other  adjustable 
parameters. 


Figure  7  inverse  correlation  length  as  a  function  of 
normalized  temperature.  Circles  are  from  neutron¬ 
scattering  measurements  of  LaMnCL 1(1  Squares  (with 
error  bars)  are  from  our  simulations  of  the  AF  Ising 
model  with  a  nonlinear  correction  to  Boltzmann 
statistics  from  3x3x3  regions.  The  3x3x3  regions  were 
chosen  to  minimize  the  energy  (as  shown  in  the  inset  of 
Fig.  8)  with  no  other  adjustable  parameters.  Conversion 
of  the  simulation  results  to  physical  units  is  made  using 
0.40  nm  for  the  average  lattice  parameter  of  LaMnCb,  as 
measured  by  neutron  scattering 
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VI.  Additional  accomplishments 

Finite-size  thermal  effects  are  found  to  significantly  influence  the  fluctuations  and 
reaction  rates  of  many  materials.  Another  example  is  the  slow  relaxation  of  disordered 
systems.  We  have  succeeded  in  simulating  several  features  in  the  dynamics  of  viscous 
liquids,  including  the  non- Arrhenius  activation  and  stretched-exponential  relaxation.  We 
find  that  the  slow  dynamics  comes  from  a  distribution  of  sharp  jumps,  as  shown  in  Fig.  8. 
The  jumps  come  from  “entropic  entrapment,”  where  the  nonlinear  correction  to 
temperature  enhances  the  stability  of  low-entropy  states,  until  the  region  fluctuates  into  a 
high-entropy  state  that  facilitates  rapid  relaxation,  somewhat  like  an  avalanche.  Direct 
evidence  for  this  picture  of  complex  relaxation  comes  from  NMR  measurements  of  jump- 
angle  distributions  (Fig.  8  inset).11  Understanding  the  dynamics  of  complex  systems  helps 
guide  a  search  for  ways  to  enhance  the  stability  of  disordered  materials  by  decreasing 
their  local  entropy.  Maintaining  low  entropy  will  help  avoid  outlier  points  exhibiting  fast 
relaxation,  which  may  act  as  “hot  spots”  that  initiate  unwanted  chemical  reactions. 

In  our  research  we  use  simplified  models  that  allow  us  to  focus  on  fundamental 
issues  in  statistical  thermodynamics.  We  seek  to  understand  the  most  basic  laws  that  yield 
universal  features  in  the  response  of  complex  systems.  We  address  the  question:  what 
minimal  rules  capture  the  essence  of  complex  dynamics?  Most  of  the  work  during  the 
previous  funding  period  involved  the  thermal  and  dynamic  properties  of  the  Ising  model, 


which  is  the  simplest  model  for  a  thermodynamic  phase  transition.  Figure  9  shows  an 
Arrhenius  plot  of  the  average  time  for  a  biased  Ising  model  to  activate  over  an  energy 


barrier,  from  a  state  of  high  potential  energy  to  the  ground  state.  Note  that  the  simplicity 
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Figure  8  Non-exponential  relaxation  as  a  function  of  time 
after  removing  an  applied  field  in  a  simulation  of  a 
supercooled  liquid.  The  solid  (black)  curve  shows  the  net 
alignment  from  81  regions,  each  of  size  16x16x16.  The 
dashed  (red)  curve  shows  a  fit  to  the  net  behavior  using  the 
stretched-exponential  function.  The  symbols  show  the 
response  of  10  individual  regions,  showing  sharp  jumps 
and  steps  in  the  local  response.  The  relatively  static  periods 
come  from  entropic  entrapment,  when  the  region  is  in  a 
state  of  low  entropy.  The  large  jumps  arise  when  the  region 
becomes  activated  by  high  entropy.  The  inset  shows  results 
from  NMR  measurements  of  the  molecular  rotation  angles 
in  a  supercooled  liquid, 1 1  showing  similar  small-angle 
diffusion  combined  with  occasional  large-angle  jumps. 
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Figure  9.  Arrhenius  plot,  showing  the  logarithm  of  the 
activation  time  as  a  function  of  the  inverse  of  the  initial 
temperature.  The  solid  lines  are  from  our  simulations  of 
the  Ising  model  (canonical  ensemble),  and  symbols  are 
from  simulations  of  the  Creutz  model  (microcanonical 
ensemble).  The  different  colors  are  for  different  system 
sizes,  as  given  in  the  legend.  The  activation  time  is 
determined  by  averaging  multiple  simulations  of  the 
time  taken  for  the  initial  (metastable)  state  to  activate 
over  an  energy  barrier  into  the  ground  state,  causing  a 
sudden  release  of  energy.  Each  data  point  comes  from 
the  average  of  several  simulations. 


of  the  Ising  model  allows  simulations  over  more  than  10  orders  of  magnitude  in  time, 
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comparable  to  the  range  of  experimental  time  scales  for  Arrhenius  activation  in  HMX. 
We  continue  to  use  Monte-Carlo  simulations  of  the  Ising  model  to  efficiently  explore  the 
minimal  complexity  needed  for  realistic  dynamics. 

Also  shown  in  Fig.  9  are  the  activation  rates  found  from  the  Creutz  model.  The 
Creutz  model  starts  with  the  interaction  energy  of  the  Ising  model,  but  adds  a  source  of 
kinetic  energy  so  that  the  simulations  can  be  done  in  the  microcanonical  ensemble.  This 
ensemble,  which  uses  conservation  of  energy  instead  of  fixed  temperature,  is  crucial  for 
studying  transient  behavior  using  a  self-consistent  local  temperature  for  systems  that  are 
out  of  equilibrium.  Figure  10  shows  an  example  of  how  the  Creutz  model  gives  the 
temperature  as  a  function  of  time  as  the  system  activates  out  of  the  high-energy  initial 
state.  The  inset  shows  details  about  the  activation  process.  Note  the  brief  reduction  in 
temperature  (kinetic  energy)  just  prior  to  the  initiation  of  activation,  as  the  system 
surmounts  the  energy  barrier  before  releasing  its  excess  potential  energy.  The  Creutz 
model  is  the  simplest  way  to  simulate  microcanonical  processes,  thus  providing  an 
elementary  link  between  Monte-Carlo  and  molecular-dynamics  simulations. 


Figure  10.  Semi-log  plot  of  temperature  as  a  function  of  time. 
The  data  points  are  from  our  simulations  of  the  Creutz  model 
using  the  microcanonical  ensemble.  Each  color  represents  a 
different  repeat  of  the  same  simulation,  starting  from  a 
statistically  different  point,  showing  the  range  of  stability  in 
the  initial  state.  The  inset  shows  details  of  the  energy  release. 
Note  the  sharpness  of  the  exothermic  activation,  and  the  brief 
reduction  in  temperature  (kinetic  energy)  that  precedes  the 
energy  release  as  the  system  surmounts  the  energy  barrier.  In 
the  inset,  the  simulations  have  been  offset  in  time  for  clarity. 
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Practical  application  of  our  fundamental  ideas  is  most  efficiently  accomplished  by 
collaboration.  Recently  I  have  collaborated  with  Professors  Don  Thompson  and  Tommy 
Sewell  at  the  University  of  Missouri,  who  are  experts  in  detailed  simulations  of  energetic 
materials  using  realistic  models.  Through  many  e-mails,  two  lengthy  teleconference  calls, 
and  four  visits  to  Columbia,  we  have  discussed  the  evidence  for  nonlinear  corrections  to 
temperature,  and  the  possible  importance  for  energetic  materials.  We  developed  a  plan  on 
how  to  proceed.  1)  Run  a  test  simulation  of  liquid  nitro methane  to  establish  the  relevance 
of  nonlinear  corrections  for  MD  simulations  of  energetic  materials,  as  shown  in  Figs.  1 
and  2.  2)  Investigate  nonlinear  corrections  to  MD  simulations  for  simplified  models, 
starting  with  the  Lennard- Jones  fluid  as  shown  in  Fig.  3,  then  transfer  the  ideas  to  more- 
sophisticated  models.  3)  Explore  models  that  bridge  the  scale  of  complexity  between  the 
Ising  model  and  models  for  energetic  materials,  such  as  the  Potts  model  for  multiple 
degrees  of  freedom  on  a  lattice.  4)  Evaluate  various  moments  in  the  distribution  of  jump 
rates,  focusing  on  outlier  statistics,  as  shown  in  Fig.  8.  5)  Maintain  our  collaboration  via 
frequent  discussions  and  timely  visits.  In  summary,  we  will  continue  to  use  computer 
simulations  to  investigate  nonlinear  corrections  to  temperature  that  are  necessary  to  fully 
describe  the  thermal  properties  of  complex  systems,  including  fluctuations  and  hot  spots 
that  may  significantly  influence  the  local  reaction  rates  in  energetic  materials. 
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